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Abstract. X-ray tomographic image reconstruction consists of determining an object 
function from its projections. In many applications such as non-destructive testing, we 
look for a fault region (air) in a homogeneous, known background (metal). The image 
reconstruction problem then becomes the determination of the shape of the default re- 
gion. Two approaches can be used: modeling the image as a binary Markov random field 
and estimating the pixels of the image, or modeling the shape of the fault and estimating 
it directly from the projections. In this work we model the fault shape by a deformable 
polygonal disc or a deformable polyhedral volume and propose a new method for directly 
estimating the coordinates of its vertices from a very limited number of its projections. 
The basic idea is not new, but in other competing methods, in general, the fault shape 
is modeled by a small number of parameters (polygonal shapes with very small number 
of vertices, snakes and deformable templates) and these parameters are estimated either 
by least squares or by maximum likelihood methods. We propose modeling the shape 
of the fault region by a polygon with a large number of vertices, allowing modeling of 
nearly any shape and estimation of its vertices' coordinates directly from the projections 
by defining the solution as the minimizer of an appropriate regularized criterion. This 
formulation can also be interpreted as a maximum a posteriori (MAP) estimate in a 
Bayesian estimation framework. To optimize this criterion we use either a simulated an- 
nealing or a special purpose deterministic algorithm based on iterated conditional modes 
(ICM). The simulated results are very encouraging, especially when the number and the 
angles of projections are very limited. 
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1. Introduction 

Tomographic image reconstruction in non destructive testing (NDT) is recent and 
consists in determining an object f(x, y) from its projects p(r, 0): 

p{r,(f)) ~ // f {x,y)S{r — X cos (f) — y sin (f)) dx dy (1) 
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In many image reconstruction applications, especially in NDT, we know that 
/(cc, y) has a constant value ci inside a region P (fault region) and another constant 
value C2 outside that region (safe or background region), e.g. metal & air. 



ci if {x, y) e P, 
C2 elsewhere 



(2) 



The image reconstruction problem becomes then the determination of the shape 
of the fault region P. In this work, without loss of generality, we assume that 
ci = 1 and C2 = and model the shape of the object by its contour. 

There has been many works in image reconstruction dealing with this problem. 
To emphasis the originality and the place of this work, we give here a summary of 
the different approaches for this problem: 

• The first approach consists in discretizing the equation (|^) to obtain: 



where, / is the discretized values of the object f{x,y) (the pixel values of the 
image), p is values of the projection data p{r,(j)), n is a vector to represent the 
modeling and measurement errors (noise) and H the discretized Radon operator. 
Then the solution is defined as the niinimizer of a compound criterion 



where A is the regularization parameter and Q and Q has to be chosen appro- 
priately to reflect our prior knowledge on the noise and on the image. This 
is the classical regularization approach of general image reconstruction problem. 
One can also interpret J(/) as the maximum a posteriori (MAP) criterion in the 
Bayesian estimation framework where exp Q(/)] represents the likelihood term 
and exp [— ri(/)] the prior probability law. 

This approach has been used with success in many applications (e.g. [|[ ||, ^, ^) 
but the cost of its calculation is huge due to the great dimension of /. Many works 
have been done on choosing appropriate regularization functionals or equivalently 
appropriate prior probability laws for / to enforce some special properties of the 
image such as smoothness, positivity or piecewise smoothness f^, 0, |[ §| . Among 
these, one can mention mainly two types of functions for ^l{f)'. 

Entropic laws: 



p = Hf + n 



(3) 



J{f)^Q{p-Hf) + xn{f), 



(4) 



N 




(5) 



Markovian laws: 



N 




(6) 



j=l igAG- 
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with convex potential functions: 



4>ix,y) = <{ \x-y\^, -|a;-y|log-, logcosh|x-y| 



or non convex potential functions: 



,y) = {min{|x-y|M},^^^^,...} 



(7) 



(8) 



See for example Q for the entropic laws, J? 
with convex potential functions and M, nC 



for scale invariant markovian laws 
[lli for markovian laws with non 



convex potential functions and other specific choices. 

• The second approach consists in modeling directly the closed contour of the 
object as the zero-crossing of a smooth function u{x, y): 



dD^{{x,y):u{x,y) = 0} and /(:z:,y)=|^^ ^^l\x,y)<0 ' 

and in defining a time evolution for u (consequently the corresponding contour dD 
and function /) such that 

dD{t)^{ix,y):u{x,y,t)^0} (10) 

and such that as time i i— s- cx) we arrive at a function u(x, y) such that the associated 
/(x, y) is a solution to the inverse problem in a least square (LS) sense. This means 
that the evolution of u and consequently the corresponding contour dD and object 
/ is such that the LS criterion 

J(/) = |b-i?(/)f (11) 

decreases during the evolution. In this approach the function u is assimilated to 
a surface (of heat front wave) and the evolution of the contours dD is constraint 
to be perpendicular to the surface. This means that the variation of the interior 
region (Sx, Sy) is such that 

Vu 

{Sx, Sy) = a(x, y, t)—— (12) 
|Vm 

This is the Level-Set approach originally developed by Osher and Sethian [|l2j for 
problems involving the motion of curves and surfaces and then adapted, used and 
referred as snakes or active contour models by many authors in computer vision 
p^ , Q and recently for inverse problems [|l5j . See also jl^, |l^ for the application 
of the deformable templates in tomographic image reconstruction. 

This approach also needs pixel or voxel representation of the image and its 
calculation cost is huge, specially in 3D imaging systems. We are presently working 
on this approach trying to extend it for minimizing a regularized criterion in place 
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of the the LS criterion ( |ll| ) and implementing it in 2D and 3D tomographic image 
reconstruction. 

• The third approach starts by giving a parametric model for the object and then 
tries to estimate these parameters using least squares (LS) or maximum likelihood 
(ML) methods. In general, in this approach one chooses a parametric model such 
as superposition of circular or elliptical homogeneous regions to be able to relate 
analytically the projections to these parameters. For example, for a superposition 
of elliptical homogeneous regions we have: 

K 

f{x,y) =^dkfk{x - ak,y - f3k) (13) 

k=l 

where = {dk,akT Pk,ak,bk,k = 1,---,K} is a vector of parameters defining 
the parametric model of the image (density values, coordinates of the centers 
and the two diameters of the ellipses). It is then easy to calculate analytically 
the projections and the relation between the data and the unknown parameters 
becomes: 

p(r, (/)) = ft.(r, 0) + n(r, (/)) (15) 

where h{r, (p; 0) has an analytic expression in 0. The LS or the ML estimate when 
the noise is assumed to be zero mean, white and Gaussian, is then given by: 

6 = argmin{||p(r, (/)) - /i(r, 0;0)f } (16) 




This approach has also been used with success in image reconstruction [g8|, |19|, 
|2l|| . But, the range of applicability of these methods is limited to the cases where 
the parametric models are actually appropriate. 

• The fourth approach which is more appropriate to our problem of shape recon- 
struction, consists in modeling directly the contour of the object by a function, 
say g{9) in a cylindrical coordinates (p, 9) such as: 

D={{x,y):x^ + y'^g^{9)}. (17) 

The next step is then to relate the projections p(r, (p) to g{9) which, in this case 
is: 

p{r,(l))^ / (5(r-pcos((/)-6'))pdpd6'. (18) 
Jo Jq 

and finally to discretize this relation to obtain: 

p = h{g)+n (19) 
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where g represents the discretized values of g{9) defining the contour of the object 
and h(g) represents the discretized version of the nonlinear operator relating 
projection data p and g. Then, one defines the solution as the argument which 
minimizes 



where ^{g) has to be chosen appropriately to reflect some regularity property of 
the object's contour. 

In this case also one can consider J{g) as the MAP criterion with Q(g) = 
\\p — h{g)\\'^ as the likelihood term and 0(g) as the prior one. 
This approach has been used in image restoration ||2^ , but it seems to be new in 
image reconstruction applications and the proposed method in this work is in this 
category. The originality of our work is to model the contour of the object by a 
piecewise linear function which means that the object is modeled as a polygonal 
region whose vertices are estimated directly from the projection data. 

2. Proposed method 

In this paper we propose to model the contour of the object (fault region) as 
a periodic piecewise linear function or equivalently to model the shape of the 
object as a polygonal region with a great number N of vertices to be able to 
approximate any shape. Then we propose to estimate directly the coordinates 
{{xj, yj),j = 1, • • • , N} of the vertices of this polygonal region from the projection 
data. 



J{g)^\\p-h{g)\\' + Xn{g), 



(20) 



y 




{x,y) e P 
ix,y)^P 



p = h{z) + n 




Figure 1: Proposed shape reconstruction modeling. 
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The idea of modeling the shape of the object as a polygonal region is not new 
and some works have been done in image reconstruction applications. See for 
example Milanfar, Karl & Willsky |2^, but in general, in these works, 

either a hypothesis of convexity of the polygonal region has been used which is 
very restrictive in real applications or the number of vertices of the polygon is 
restricted to a very small number. In our work we do not make neither of these 
hypothesis; the polygon can be convex or not and we choose N sufficiently great 
to to be able to approximate appropriately any shape. 

The solution is then defined as the minimizer the following criterion 



J{z) = \\p - h{z)\\^ + XQ{z) 



(21) 



where z ~ x + iy is a complex vector whose real and imaginary parts represent 
the X and the y coordinates of the polygon vertices, h{z) represents the direct 
operator which calculates the projections for any given z and ^l{z) is chosen to be 
a function which reflects the regularity of the object contour. One can for example 
choose the following function: 



N 



(22) 



which favors a shape whose contour length is near a prior known value 27ri?o- In 
this work we used the following function: 



N 



N 



Zj+1 



)/2| 



(23) 



which favors a shape whose local curvature is limited. Note that \zj — {zj-i — Zj+i)/2\ 
is just the Euclidian distance between the point Zj and the midpoint of the line seg- 
ment passing through Zj-i and Zj+i and so this choice favors a shape whose local 
curvature is limited. We can also give a probabilistic interpretation to this choice. 
In fact we can consider Zj as random variables with the following Markovian law: 



p{zj\z) = p{zj\zj^i,Zj+i) cx exp 



(24) 



Other functions are possible and are studied in this work. 

In both cases, the criterion J{z) is multi-modal essentially due to the fact that 
h{z) is a nonlinear function of z. Calculating the optimal solution corresponding 
to the global minimum of ( ^ ) needs then carefully designed algorithms. For this 
we propose the two following strategies: 

• The first is to use a global optimization technique such as simulated annealing. 
This technique has given satisfactory results as it can be seen from the simulations 
in the next section. However, this algorithm needs a great number of iterations 
and some skills for choosing the first temperature and cooling schedule, but the 
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overall calculations is not very important due to the fact that, in this algorithm, 
at each iteration only one of the vertices Zj is changed. So, at each step we need 
to calculate the variation of the criterion due to this change which can be done 
with a reasonable cost. 

• The second is to find an initial solution in the attractive region of the global 
optimum and to use a local descent type algorithm such as the ICM (Iterated 
conditional modes) of Besag |2^, |2^, |2|] to find the solution. 

The main problem here is how to find this initial solution. For this, we used 
a moment based method proposed by Milanfar, Karl & Willsky ^ ^ which 
is accurate enough to obtain an initial solution which is not very far from the 
optimum. The basic idea of this method is to relate the moments of the projections 
to the moments of a class of polygonal regions obtained by an affine transformation 
of a regular polygonal region, and so to estimate a polygonal region whose corners 
are on an ellipse and whose moments up to the second order matches those of the 
projections. 

However, there is no theoretical proof that this initial solution will be in the 
attractive region of the global optimum. The next section will show some results 
comparing the performances of these two methods as well as a comparison with 
some other classical methods. 



3. Simulation results 

To measure the performances of the proposed method and keeping the objective 
of using this method for NDT applications where the number of projections are 
very limited, we simulated two cases where the objects have polygonal shapes with 
= 40 corners (hand-made) and calculated their projections for only 5 directions 

= {-45,-22.5,0,22.5,45} degrees. 

Then, we added some noise (white, Gaussian and centered) on them to simulate 
the measurement errors. The signal to noise ratio (SNR) was chosen 20dB. We 
define SNR as follows: ^ 

SNR =lQ\og ^'f' 

where pi and iii are respectively data and noise samples and p and n their respective 
mean values. 

Finally, from these data we estimated the solutions by the proposed method 
using either the simulated annealing (SA) or the iterated conditional modes (ICM) 
algorithms. 

Figure 2 shows these two objects and their relative simulated projections data. 

In Figures 3 and 4, we give the reconstruction results obtained by the proposed 
method using either the S A algorithm (Figure 3) or the ICM algorithm (Figure 4) . 
In these figures we show the original objects, the initial solutions, the intermediate 
solutions during the iterations and the final reconstructed objects obtained after 
200 iterations. 
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Note that, the SA algorithm is theoretically independent of initialization while 
the ICM is not. However, in these figures, for the purpose of the comparison, we 
show the results obtained by the two algorithms with the same initialization. 

To show that the method is not very sensible on the prior knowledge of the 
vertices number, we give in Figure 5, the reconstruction results of the object 2 in 
4 cases: = 10, iV = 20 and = 30 and A^ = 40. As we can remark all the 
reconstructed results seem satisfactory. 

In Figure 6 we show a comparison between the results obtained by the proposed 
method and those obtained either by a classical backprojection method or by some 
other methods in the first approach using (H) and (jij) with different regularization 
functionals ^{f), more specifically: 

— Gaussian Markov models (|^) with the potential function (j){x,y) = |a; — 
which can also be considered as a quadratic regularization method; and 

— Compound Markov models with non convex potential functions 4'{x, y) = 
min {|a; — j/p, l} which is a truncated quadratic potential function. 

In the first case the criterion to optimize is convex and we used a conjugate gradient 
(CG) algorithm to find the optimized solution. In the second case the criterion 
is not convex and we used a Graduated non convexity (GNC) based optimization 
algorithm developed in jl^, ||, ^ to find the solution. 

Note that, these results are given here to show the relative performances of 
these methods in a very difficult situation where we have only five projections. 
In fact, in more comfortable situations (more projections uniformly distributed 
around the object and high SNR) all these methods, even the simplest one such 
as the classical backprojection will give similar and satisfactory results. Here, we 
compare the results obtained from the same restricted set of data. 



8 



Projections 





Figure 4: Reconstruction results using a moment-based initialization and a local 
descent (ICM) minimizer. 

o) Original objects, +) Initializations, .) Evolution of the solutions during 

the iterations and 

•k) Final reconstructed objects. 
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Figure 6: A comparison with backprojection and some other classical methods: 

a) Original objects, 

b) Results obtained by the proposed method using the SA optimization algorithm, 

c) Results obtained by the proposed method using the ICM algorithm, 

d) Backprojection, 

e) Gaussian Markov modeling MAP estimation and 

f) Compound Markov modeling and GNC optimization algorithm. 
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Figure 7: A comparison with backprojection and some other classical methods: 

a) Original objects, 

b) Results obtained by the proposed method using the SA optimization algorithm, 

c) Results obtained by the proposed method using the ICM algorithm, 

d) Backprojection, 

e) Gaussian Markov modeling MAP estimation and 

f) Compound Markov modeling and GNC optimization algorithm. 
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4. Conclusions 

A new method for tomographic image reconstruction of a compact binary object 
from a small number of its projections is proposed. The basic idea of the proposed 
method is to model the compact binary object as a polygonal region whose vertices 
coordinates are estimated directly from the projections using the Bayesian MAP 
estimation framework or equivalently by optimizing a regularized criterion. 

Unfortunately, this criterion is not unimodal. To find the optimized solution 
two algorithms have been proposed: 

- a global optimization algorithm based on simulated annealing (SA) and 

- a local descent-based method based on the Iterated Conditional Modes (ICM) 
algorithm proposed originally by Besag, with a good initialization obtained by 
using a moment based method. 

The first algorithm seems to give entire satisfaction. The second can also give 
satisfaction, but it may also be plugged in a local minimum. In both algorithms 
the main cost calculation is due to the calculus of the variation of the criterion 
when one of the vertices coordinates is changed. We have written an efhcient 



program to do this 129, 30 



An extension of this work to 3D image reconstruction with small number of 
conic projections is in preparation | |3l| , The final objective of the proposed 
method is for non destructive testing (NDT) image reconstruction applications 
where we can use not only X-rays but also ultrasound or Eddy currents or a 
combination of them ||3^, Q |l^ to localize and to characterize more accurately 
any anomalies (air bulbs) in metallic structures. 
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